In [1]:

    

import numpy as np

import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 8.0)    # resize plots

from scipy import stats

import blueice as bl
from blueice.test_helpers import test_conf


    

    




C:\Anaconda3\lib\site-packages\matplotlib\__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)

In [2]:

    

m = bl.Model(test_conf())
m.expected_events()


    

    
Out[2]:





array([1000])

In [3]:

    

q = np.linspace(-5, 5, 100)
plt.plot(q, m.sources[0].pdf(q))


    

    
Out[3]:





[<matplotlib.lines.Line2D at 0xa3b5358>]






    

In [4]:

    

lf = bl.UnbinnedLogLikelihood(test_conf())
lf.add_shape_parameter('some_multiplier', (0.5, 1, 2))
lf.add_rate_parameter('s0')
lf.prepare()
d = lf.base_model.simulate()
print(len(d), lf.base_model.expected_events())
lf.set_data(d)


    

    










    




987 [1000]

In [5]:

    

for k, m in lf.anchor_models.items():
    print(k, m.expected_events(), m.sources[0].events_per_day, m.sources[0].fraction_in_range)


    

    




(0.5,) [ 500.] 500.0 1
(1.0,) [1000] 1000 1
(2.0,) [2000] 2000 1

In [6]:

    

x = np.linspace(0.5, 2, 100)
lf.plot_likelihood_ratio(('s0_rate_multiplier', x), ('some_multiplier', x), vmax=100)
plt.xlim(0.5, 2)
plt.ylim(0.5, 2)


    

    










    
Out[6]:





(0.5, 2)






    

In [7]:

    

# Should be flat 0 (since fitting foo_rate_multiplier completely compensates for some_multiplier)
lf.plot_likelihood_ratio(('some_multiplier', x))

# Should be centered at 1 (approximately, depends on len(d)):
lf.plot_likelihood_ratio(('some_multiplier', x), s0_rate_multiplier=1)

# Should be centered at 0.5:
lf.plot_likelihood_ratio(('some_multiplier', x), s0_rate_multiplier=2)

# Should be centered at 2
lf.plot_likelihood_ratio(('some_multiplier', x), s0_rate_multiplier=0.5)
plt.ylim(-1, 5)


    

    
Out[7]:





(-1, 5)






    

In [8]:

    

# Should be flat 0....
# Except when len(d) deviates significantly from 1000, 
# in which case you'll see it following the light blue/red curve.
# Unlike the rate multiplier, the shape parameter is bounded by the anchors we put in, so it can't take a value
# outside (0.5, 2).
lf.plot_likelihood_ratio(('s0_rate_multiplier', x))

# Should be centered at 1
lf.plot_likelihood_ratio(('s0_rate_multiplier', x), some_multiplier=1)

# Should be centered at 0.5
lf.plot_likelihood_ratio(('s0_rate_multiplier', x), some_multiplier=2)

# Should be centered at 2:
lf.plot_likelihood_ratio(('s0_rate_multiplier', x), some_multiplier=0.5)
plt.xlim(0.5, 2)
plt.ylim(-1, 5)


    

    
Out[8]:





(-1, 5)






    

In [9]:

    

# Will be centered at 2, and very steep. 
# It would like to go to 10, but can't because of the bounds on the shape param
lf.plot_likelihood_ratio(('some_multiplier', np.linspace(0, 3, 100)), s0_rate_multiplier=0.1)


    

In [10]:

    

# Should be exactly the same:
plt.plot(x, [lf(some_multiplier=q) for q in x])
plt.plot(x, [lf(s0_rate_multiplier=q) for q in x])


    

    
Out[10]:





[<matplotlib.lines.Line2D at 0xab26780>]






    

In [11]:

    

w1 = lf.make_objective(some_multiplier=0.6)
w2 = lf.make_objective(s0_rate_multiplier=0.6)
print(w1, "\n", w2)
f1 = w1[0]
f2 = w2[0]


    

    




(<function make_objective.<locals>.objective at 0x00000000090CE488>, ['s0_rate_multiplier'], array([1]), [(0, None)]) 
 (<function make_objective.<locals>.objective at 0x000000000A574D90>, ['some_multiplier'], array([1]), [(0.5, 2)])

In [12]:

    

# Should be exactly the same
plt.plot(x, [f1([q]) for q in x])
plt.plot(x, [f2([q]) for q in x])


    

    
Out[12]:





[<matplotlib.lines.Line2D at 0xa601eb8>]






    

In [13]:

    

q = 0.6
kwargs = dict()
a = lf.bestfit_minuit(some_multiplier=q, 
                      minimize_kwargs=dict(print_level=0),
                      **kwargs)[0]
b = lf.bestfit_minuit(s0_rate_multiplier=q, **kwargs)[0]
print(a, b, len(d)/lf.base_model.expected_events()/q)


    

    




{'s0_rate_multiplier': 1.6449356394835803} {'some_multiplier': 1.6449356394854142} [ 1.645]